import geopandas as gpd
from openpyxl import Workbook, load_workbook
from tqdm import tqdm
import fiona
import pandas as pd
import requests
from shapely.geometry import MultiPolygon,Polygon, MultiLineString, LineString
from pyproj import CRS
import matplotlib.pyplot as plt
from bs4 import BeautifulSoup
import concurrent.futures
from urllib.parse import urljoin, urlparse


def shp2xlsData(path, bdf_path):
    # with fiona.open(path) as src:
    #     print("SHP文件完整架构:", src.schema)
    #     print("所有列名:", list(src.schema['properties'].keys()))

    # bdf_data = pd.read_csv(bdf_path,encoding='gb2312')

    # bdf_data.to_excel("PSR_84_New.xlsx",
    #                   index=False,
    #                   engine='openpyxl'
    #                   )

    # print("属性表数据成功写入。")

    wb = Workbook()
    wv = wb.active

    wv.append(['gid', 'area_km2', 'name', 'lon', 'lat', "Upper Left Longitude",
               "Upper Left Latitude", "Upper Right Longitude", "Upper Right Latitude", "Lower Right Longitude",
               "Lower Right Latitude",
               "Lower Left Longitude", "Lower Left Latitude"
               ])

    shapefile = gpd.read_file(path)
    shapefile.info()

    # 遍历所有的行列
    # rows, columns = shapefile.shape

    selected_columns = shapefile.columns.tolist()
    print(selected_columns)
    # selected_columns = ['gid', 'area_km2', 'name', 'CENTROID_X', 'CENTROID_Y']
    for _, row in tqdm(shapefile.iterrows(), total=len(shapefile), desc="写入Excel"):
        row_data = [row[col] for col in selected_columns]

        # print(row_data[5].bounds)
        # 组一下最大外包矩形框的坐标
        minx, miny, maxx, maxy = row_data[5].bounds

        point = [
            (minx, miny),  # 左下角
            (maxx, miny),  # 右下角
            (maxx, maxy),  # 右上角
            (minx, maxy)  # 左上角
        ]

        for tuple_tmp in point:
            row_data.append(tuple_tmp[0])
            row_data.append(tuple_tmp[1])

        data = [value for index, value in enumerate(row_data) if index != 5]
        wv.append(data)

    wb.save("PSR_84_New.xlsx")
    print("已经将shp数据转换为xlsx")


def getShadowId(shadow_id_path):
    # 先读取shadow文件里面的数据
    wb = load_workbook(shadow_id_path, read_only=True)

    ws = wb.active
    shadow_id_array = []
    flag = True
    for row in ws.rows:
        if flag:
            flag = False
            continue
        row_data = []
        for _, cell in enumerate(row):
            row_data.append(cell.value)
        shadow_id_array.append(row_data)

    return shadow_id_array


def PutShadowID(shadow_id_array, xlsx_path, Numbers):
    # print(shadow_id_array)

    # 怕网页标签元素
    base_url = "https://data.im-ldi.com/mds/shadowcam_published/"
    headers = {
        "User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/91.0.4472.124 Safari/537.36",
        "Accept": "text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,*/*;q=0.8",
        "Accept-Language": "zh-CN,zh;q=0.9,en;q=0.8",
        "Connection": "keep-alive",
    }

    wb = Workbook(
    )
    wv = wb.active

    wv.append(["影像编号",
               "Center Latitude",
               "Center Longitude",
               "Lower Left Latitude",
               "Lower Left Longitude",
               "Lower Right Latitude",
               "Lower Right Longitude",
               "Upper Left Latitude",
               "Upper Left Longitude",
               "Upper Right Latitude",
               "Upper Right Longitude"
               ])

    for idx, shadow in enumerate(tqdm(shadow_id_array, total=(Numbers[1] - Numbers[0] + 1), desc="爬取影像进度：")):

        if idx not in range(Numbers[0], Numbers[1] + 1):
            continue
        # base_url+id为目的地址
        url = base_url + shadow[0]

        try:
            response = requests.get(url, headers=headers)
            if response.status_code != 200:
                print("无法访问！")
                return

            # 爬取元素
            # 解析HTML内容
            soup = BeautifulSoup(response.text, 'html.parser')
            count = 0
            localtion_temp = [shadow[0]]

            for tbody in soup.find_all('tbody'):
                # print(tbody)
                # print(count)
                # 第7个就是需要的table
                if count == 7:
                    td_data = [td.get_text(strip=True) for td in tbody.find_all('td')]
                    print(td_data)

                    for index, value in enumerate(td_data):

                        if index == 10:
                            break
                        localtion_temp.append(float(value))
                        # localtion_dict[Text_array[index]] = float(value)
                    # print(localtion_dict)

                count += 1
        except Exception as e:
            print("爬取网页时候出现错误！")
            return

        wv.append(localtion_temp)

    wv.save("getURL_Location_Data.xlsx")
    print("网页数据爬取完成！")

    return


def process_shapefile(shp_path, location_data_path, url_data_path, prj_path, lat_min=-81, area_over=0.2):
    # 读取shapefile
    shapefile = gpd.read_file(shp_path)

    # 读取Location数据
    location_data = pd.read_excel(location_data_path)

    # 读取URL数据
    url_data = pd.read_excel(url_data_path)

    # 创建一个新的DataFrame用于保存结果
    output_data = []

    # 读取moon_2000_slope.prj文件并获取投影信息
    with open(prj_path, 'r') as prj_file:
        prj_text = prj_file.read().strip()
    crs_target = CRS.from_string(prj_text)  # 读取目标投影
    # 获取shapefile原始坐标系
    original_crs = shapefile.crs

    # 将shapefile转换为目标坐标系（moon_2000_slope）
    shapefile = shapefile.to_crs(crs_target)  # 转换为目标坐标系

    # 遍历shapefile的每一行
    for idx, row in shapefile.iterrows():
        # 获取shapefile的几何边界
        geom = row['geometry']
        # 获取shapefile的gid、area_km2和name以及质心坐标
        gid = row['gid']
        area_km2 = row['area_km2']  # 转换为平方千米
        name = row['name']
        centroid_x, centroid_y = geom.centroid.x, geom.centroid.y
        # 获取边界（即外部的LineString）
        boundary = geom.boundary

        # 如果边界是一个LineString（有时会是MultiLineString），则可以将它转换为一个多边形
        if isinstance(boundary, MultiLineString):
            # 选择最大的线段（或根据其他标准选择），例如最大长度的线段
            largest_line = max(boundary.geoms, key=lambda line: line.length)
            # 创建一个闭合的Polygon（假设该线段能够封闭成一个多边形）
            coords = list(largest_line.coords)
            coords.append(coords[0])  # 闭合线段
            geom = Polygon(coords)  # 将其转换为Polygon

        elif isinstance(boundary, Polygon):
            # 如果直接是Polygon类型，直接处理
            geom = boundary

        # 遍历Location_Data_path中的每行，检查是否与shapefile相交
        for _, loc_row in location_data.iterrows():
            # 读取Location的四个角坐标
            loc_center_lat = loc_row['Center Latitude']
            loc_lower_left_lon = loc_row['Lower Left Longitude']
            loc_lower_left_lat = loc_row['Lower Left Latitude']
            loc_lower_right_lon = loc_row['Lower Right Longitude']
            loc_lower_right_lat = loc_row['Lower Right Latitude']
            loc_upper_left_lon = loc_row['Upper Left Longitude']
            loc_upper_left_lat = loc_row['Upper Left Latitude']
            loc_upper_right_lon = loc_row['Upper Right Longitude']
            loc_upper_right_lat = loc_row['Upper Right Latitude']
            if loc_center_lat >= lat_min:
                continue  # 如果纬度不符合条件，跳过当前记录
            # 将location的四个交点转换到目标坐标系
            loc_polygon = Polygon([
                (loc_upper_left_lon, loc_upper_left_lat),
                (loc_upper_right_lon, loc_upper_right_lat),
                (loc_lower_right_lon, loc_lower_right_lat),
                (loc_lower_left_lon, loc_lower_left_lat),
            ])

            # 将location polygon转换为与shapefile相同的坐标系
            loc_polygon_gdf = gpd.GeoDataFrame(geometry=[loc_polygon], crs=original_crs)  # 使用shapefile的原始坐标系
            loc_polygon_gdf = loc_polygon_gdf.to_crs(crs_target)  # 转换为目标投影

            # 判断shapefile的几何与Location四边形是否相交
            intersection_found = geom.intersects(loc_polygon_gdf.geometry[0])

            if intersection_found:
                # 查找对应影像编号的URL信息
                image_id = loc_row['影像编号']
                url_row = url_data[url_data['影像编号'] == image_id]
                if not url_row.empty:
                    image_url = url_row.iloc[0]['URL']  # 获取URL列数据
                    # 如果有相交，记录结果
                    output_data.append({
                        'gid': gid,
                        'area_km2': area_km2,
                        'name': name,
                        'CENTROID_X': centroid_x,
                        'CENTROID_Y': centroid_y,
                        '影像编号': loc_row['影像编号'],
                        'URL': image_url  # 保存URL
                    })

    return output_data


def process_location(loc_row, geom, lat_min, crs_target, original_crs, url_data, gid, area_km2, name, centroid_x,
                     centroid_y):
    # 读取Location的四个角坐标
    loc_center_lat = loc_row['Center Latitude']
    loc_lower_left_lon = loc_row['Lower Left Longitude']
    loc_lower_left_lat = loc_row['Lower Left Latitude']
    loc_lower_right_lon = loc_row['Lower Right Longitude']
    loc_lower_right_lat = loc_row['Lower Right Latitude']
    loc_upper_left_lon = loc_row['Upper Left Longitude']
    loc_upper_left_lat = loc_row['Upper Left Latitude']
    loc_upper_right_lon = loc_row['Upper Right Longitude']
    loc_upper_right_lat = loc_row['Upper Right Latitude']

    if loc_center_lat >= lat_min:
        return None  # 如果纬度不符合条件，跳过当前记录

    # 将location的四个交点转换到目标坐标系
    loc_polygon = Polygon([
        (loc_upper_left_lon, loc_upper_left_lat),
        (loc_upper_right_lon, loc_upper_right_lat),
        (loc_lower_right_lon, loc_lower_right_lat),
        (loc_lower_left_lon, loc_lower_left_lat),
    ])

    # 将location polygon转换为与shapefile相同的坐标系
    loc_polygon_gdf = gpd.GeoDataFrame(geometry=[loc_polygon], crs=original_crs)  # 使用shapefile的原始坐标系
    loc_polygon_gdf = loc_polygon_gdf.to_crs(crs_target)  # 转换为目标投影

    # 判断shapefile的几何与Location四边形是否相交
    intersection_found = geom.intersects(loc_polygon_gdf.geometry[0])

    if intersection_found:
        # 查找对应影像编号的URL信息
        image_id = loc_row['影像编号']
        url_row = url_data[url_data['影像编号'] == image_id]
        if not url_row.empty:
            image_url = url_row.iloc[0]['URL']  # 获取URL列数据
            # 如果有相交，记录结果
            return loc_polygon_gdf.geometry[0], image_id, image_url  # 返回影像编号及其几何体
    return None  # 没有相交的返回None


def select_minimal_covering(idx,geom, images, output_prefix="coverage_step"):
    """
    贪心选择最少影像覆盖目标几何，每选择一张影像时保存一张所有已选影像交集的可视化图。
    """
    # 按交集面积降序排列
    images = sorted(images, key=lambda x: x[0].intersection(geom).area, reverse=True)

    covered_area = Polygon()
    selected_images = []

    target_gdf = gpd.GeoSeries([geom])

    # 先选最大面积影像
    if not images:
        print("No intersecting images found. Check geom or image geometries.")
        return []
    largest_geom, largest_id, largest_url = images[0]
    largest_intersection = largest_geom.intersection(geom)
    covered_area = covered_area.union(largest_intersection)
    selected_images.append((largest_id, largest_url, largest_intersection))

    # 可视化当前所有 selected_images 的覆盖
    # _save_selected_images_plot(target_gdf, selected_images, step=largest_id, prefix=output_prefix)

    i = 1
    while i < len(images):
        current_geom, current_id, current_url = images[i]
        current_intersection = current_geom.intersection(geom)

        if covered_area.contains(current_intersection.buffer(-1e-8)):
            i += 1
            continue
        uncovered_part = current_intersection.difference(covered_area)
        replaced = False
        for j in range(i + 1, len(images)):
            next_geom, next_id, next_url = images[j]
            next_intersection = next_geom.intersection(geom)
            # _visualize_checking_replacement(target_gdf, uncovered_part, next_intersection)
            if next_intersection.contains(uncovered_part.buffer(-1e-8)):
                covered_area = covered_area.union(next_intersection)
                selected_images.append((next_id, next_url, next_intersection))
                # _save_selected_images_plot(target_gdf, selected_images, step=next_id, prefix=output_prefix)
                replaced = True
                break

        if not replaced:
            covered_area = covered_area.union(current_intersection)
            selected_images.append((current_id, current_url, current_intersection))
            # _save_selected_images_plot(target_gdf, selected_images, step=current_id, prefix=output_prefix)

        if covered_area.contains(geom.buffer(-1e-8)):
            break

        i += 1

    print("FID: ", idx)
    print("Final uncovered area (geom - covered_area): ", geom.difference(covered_area).area)
    return selected_images


def _visualize_checking_replacement(target_gdf, uncovered_part, next_intersection, step=1):
    """
    可视化正在判断的 uncovered_part 和候选 next_intersection
    """
    fig, ax = plt.subplots(figsize=(8, 8))

    # 目标区域
    target_gdf.plot(ax=ax, facecolor='none', edgecolor='blue', linewidth=1.5, label='Target Geometry')

    # 未覆盖部分（当前影像中尚未被覆盖）
    gpd.GeoSeries([uncovered_part]).plot(ax=ax, color='red', alpha=0.4, edgecolor='darkred', label='Uncovered Part')

    # 替代候选影像的交集区域
    gpd.GeoSeries([next_intersection]).plot(ax=ax, color='green', alpha=0.4, edgecolor='darkgreen',
                                            label='Next Intersection')

    ax.set_title(f"Check if Step {step} Can Cover Uncovered Part")
    ax.legend()
    plt.axis('equal')
    plt.tight_layout()
    plt.savefig(f"check_replace_step_{step}.png", dpi=150)
    plt.close(fig)


def _save_selected_images_plot(target_gdf, selected_images, step, prefix):
    """
    可视化当前所有 selected_images 中的交集区域，并保存图像。
    """
    fig, ax = plt.subplots(figsize=(8, 8))

    # 绘制目标区域
    target_gdf.plot(ax=ax, facecolor='none', edgecolor='blue', linewidth=1.5, label='Target Geometry')

    # 绘制所有 selected_images 的交集区域
    for idx, (image_id, image_url, intersection_geom) in enumerate(selected_images):
        gpd.GeoSeries([intersection_geom]).plot(
            ax=ax,
            alpha=0.3,
            edgecolor='black',
            label=f'Image {image_id}' if idx == 0 else None  # 只让第一个图例显示一次
        )

    ax.set_title(f"Step {step}: Selected Image Coverage")
    ax.legend()
    plt.axis("equal")
    plt.tight_layout()
    plt.savefig(f"{prefix}_step_{step}.png", dpi=150)
    plt.close(fig)


def get_intersecting_shapefile(shp_path, location_data_path, url_data_path, prj_path, lat_min=-81, area_over=0.2):
    # 读取shapefile
    shapefile = gpd.read_file(shp_path)

    # 读取URL数据
    url_data = pd.read_excel(url_data_path)

    # 创建一个新的DataFrame用于保存结果
    output_data = []

    # 读取moon_2000_slope.prj文件并获取投影信息
    with open(prj_path, 'r') as prj_file:
        prj_text = prj_file.read().strip()
    crs_target = CRS.from_string(prj_text)  # 读取目标投影
    # 获取shapefile原始坐标系
    original_crs = shapefile.crs

    # # 将shapefile转换为目标坐标系（moon_2000_slope）
    shapefile = shapefile.to_crs(crs_target)  # 转换为目标坐标系
    shapefile["geometry"] = shapefile["geometry"].apply(
        lambda geom: MultiPolygon([geom]) if isinstance(geom, Polygon) else geom
    )
    # # 可视化
    # shapefile.plot()
    # plt.title("Shapefile Preview")
    # plt.axis("equal")  # 保持比例尺
    # plt.savefig(f"SHP84.png", dpi=150)
    # 读取Location数据
    location_data = pd.read_excel(location_data_path)

    # 遍历shapefile的每一行
    for idx, row in shapefile.iterrows():
        # 获取shapefile的几何边界
        geom = row['geometry']
        # 获取shapefile的gid、area_km2和name以及质心坐标
        gid = row['gid']
        area_km2 = row['area_km2']  # 转换为平方千米
        name = row['name']
        centroid_x, centroid_y = geom.centroid.x, geom.centroid.y
        # 获取边界（即外部的LineString）
        boundary = geom.boundary

        # 如果边界是一个LineString（有时会是MultiLineString），则可以将它转换为一个多边形
        if isinstance(boundary, (MultiLineString, LineString)):
            # 选择最大的线段（或根据其他标准选择），例如最大长度的线段
            largest_line = max(boundary.geoms, key=lambda line: line.length)
            # 创建一个闭合的Polygon（假设该线段能够封闭成一个多边形）
            coords = list(largest_line.coords)
            coords.append(coords[0])  # 闭合线段
            geom = Polygon(coords)  # 将其转换为Polygon

            # 使用ThreadPoolExecutor并行处理每一行location数据
            with concurrent.futures.ThreadPoolExecutor() as executor:
                futures = []
                for _, loc_row in location_data.iterrows():
                    futures.append(
                        executor.submit(process_location, loc_row, geom, lat_min, crs_target, original_crs, url_data,
                                        gid, area_km2, name, centroid_x, centroid_y))

                # 获取并处理所有结果
                intersecting_images = []
                for future in concurrent.futures.as_completed(futures):
                    result = future.result()
                    if result:
                        intersecting_images.append(result)

                # 选择最小的影像集覆盖目标几何
                selected_images = select_minimal_covering(idx,geom, intersecting_images)

                # 记录结果
                for image_id, image_url, geom in selected_images:
                    output_data.append({
                        'gid': gid,
                        'area_km2': area_km2,
                        'name': name,
                        'CENTROID_X': centroid_x,
                        'CENTROID_Y': centroid_y,
                        '影像编号': image_id,
                        'URL': image_url
                    })

                # 可视化结果
                fig, ax = plt.subplots(figsize=(10, 10))

                # 绘制shapefile的几何
                shapefile.plot(ax=ax, color='lightblue', edgecolor='black', alpha=1)

                # 绘制选定的影像
                for image_geom in [image_geom for _, _, image_geom in selected_images]:
                    gpd.GeoSeries(image_geom).plot(ax=ax, color='green', edgecolor='black', alpha=0.2)

                # 设置标题和标签
                ax.set_title(f"Shapefile {gid} and Selected Images")
                ax.set_xlabel("Longitude")
                ax.set_ylabel("Latitude")

                # 保存为PNG文件
                plt.savefig(f"shapefile_{gid}_selected_images.png", format="png")
                plt.close(fig)

    # 将结果保存为Excel
    output_df = pd.DataFrame(output_data)
    output_df.to_excel('output_result1.xlsx', index=False)

    print("处理完成，结果已保存到output_result.xlsx")


if __name__ == "__main__":
    # 转换shp格式数据为xls
    # shp 的文件位置
    # 首先要去读取PR84的shp文件找到对应坑的位置
    shp_path = "./PSR84/PSR_84.shp"
    bdf_path = "./PSR84/PSR_84.dbf"
    # shp2xlsData(shp_path,bdf_path)

    # 然后根据坑的经纬度信息去找有tif影像能覆盖到他的影像ID

    xlsx_path = "./PSR_84_New1.xlsx"
    shadow_id_path = "./shadow_id_URL_tif.xlsx"
    #
    # shadow_id_array = getShadowId(shadow_id_path)
    # # 你要爬取的影像 序号，0到666就是代表爬取前 666条数据 已知：总共有10447条数据
    # Numbers = (0, 666)
    #
    # PutShadowID(shadow_id_array, xlsx_path, Numbers)
    Location_Data_path = "./getURL_Location_Data_South.xlsx"
    prj_path = "./Moon2000_spole.prj"
    # 根据影像ID 从爬取的所有URL对应表中查询其对应的URL 放到新PSR_84.xls中
    get_intersecting_shapefile(shp_path, Location_Data_path, shadow_id_path, prj_path, -81)
